private :
      
            
      // Used to transmit submatrices between processes
      template< typename T, int Rank, bool Sparsity >
      class SubmatrixPacket;
      friend class SubmatrixPacket;
      
      // position of first entry relative to global matrix, row first; changes
      std::pair< _Index, _Index > global_starts_;
      
      // Block row and column numbers in global matrix, row first; changes
      std::pair< _Index, _Index > block_;
      
      // submatrix id used by Muesli
      int submatrixId;
      
      // host processor id set by Muesli
      int procId;
      
      // "Zero" value for data type
      _Scalar zero_;
      
       
public : 


      inline setSubmatrixId(int newid) const {
      
            submatrixId = newid;
            
      }


      inline setProcId(int newid) const {
      
            procId = newid;
            
      }


      inline const int & getSubmatrixId() const {
      
            return submatrixId;
            
      }


      inline const int & getProcId() const {
      
            return procId;
            
      }


      // Public member functions expected in MSL::Submatrix
      
      inline void debug() const {
      
            this->print();
      
      }
      
      
      inline SparseMatrix< _Scalar, _Options, _Index > * clone() const {
      
            SparseMatrix< _Scalar, _Options, _Index > * copy = new SparseMatrix< _Scalar, _Options, _Index >(*this);
            return copy;
      
      }
      
      
      inline SparseMatrix< _Scalar, _Options, _Index > * cloneForHold(int cloneId, bool swapData) const {
      
            SparseMatrix< _Scalar, _Options, _Index > * copy = new SparseMatrix< _Scalar, _Options, _Index >();
            copy->global_starts_ = this->global_starts_;
            copy->block_ = this->block_;
            if (cloneId == -1) {
                  copy->submatrixId = this->submatrixId;
            } else {
                  copy->submatrixId = cloneId;
            }
            copy->procId = this->procId;
            copy->zero_ = this->zero_;
            copy->Rank = this->Rank;
            if (swapData) {
                  this->swap(copy);
            }
            return copy;
      
      }


      inline int getColumnIndexGlobal(int k) const {
      
            return this->getColumnIndexLocal(k) + this->global_starts_.second;
      
      }
      
      
      inline int getColumnIndexLocal(int k) const {
      
            if (_Options == RowMajor) {
                  return (int) *(this->InnerIndexPtr() + k);
            } else {
                  for ((_Index) c = 0; c < this->cols(); ++c) {
                        if (*(this->OuterIndexPtr() + c + 1) > k) {
                              return (int) c;
                        }
                  }
            }
      
      }
      
      inline int getColumnIndexStart() const {
      
            return this->global_starts_.second;
      
      }
      
      
      inline bool columnIsLocal(int globalColumnIndex) const {
      
            return (this->global_starts_.second <= globalColumnIndex && globalColumnIndex < this->global_starts_.second + this->cols());
      
      }
      
      
      inline _Scalar getElement(int localRowIndex, int localColumnIndex) const {
      
            return this->coeff((_Index) localRowIndex, (_Index) localColumnIndex);
      
      }
      
      
      inline _Scalar getElementLocal(int k) const {
      
            _Index r = this->getRowIndexLocal(k);
            _Index c = this->getColumnIndexLocal(k);
            return getElement(r, c);
      
      }
      
      
      inline int getElementCount(bool includeZeros) const {
      
            if (includeZeros) {
                  return this->size();
            } else {
                  return this->nonZeros();
      
      }
      
      
      inline int getElementCount() const {
      
            return this->getElementCount(true);
      
      }


      inline int getElementCountLocal() const {
      
            return this->getElementCount(false)();
      
      }
      
      
      inline int getLocalM() const {
      
            return this->cols();
      
      }
      
      
      inline int getLocalN() const {
      
            return this->rows();
      
      }
      
      
      inline int getRowIndexGlobal(int k) const {
      
            return this->getRowIndexLocal(k) + this->global_starts_.first;
      
      }
      
      
      inline int getRowIndexLocal(int k) const {
      
            if (!(_Options == RowMajor)) {
                  return (int) *(this->InnerIndexPtr() + k);
            } else {
                  for ((_Index) r = 0; r < this->rows(); ++r) {
                        if (*(this->OuterIndexPtr() + r + 1) > k) {
                              return (int) r;
                        }
                  }
            }

      }      
      
      inline int getRowIndexStart() const {
      
            return this->global_starts_.first;
      
      }
      
      
      inline _Scalar ** getUncompressed() const {
      
            _Scalar ** result = new _Scalar [this->rows()][this->cols()];
            for (int r = 0; r < this->rows(); ++r) {
                  for (int c = 0; c < this->cols(); ++c) {
                        result[r][c] = this->coeff((_Index) r, (_Index) c);
                  }
            }
            return result;
      
      }
      
      
      inline bool isEmpty() const {
      
            return (this->nonZeros() == 0);
      
      }
      
      
      inline void pack() {
      
            this->makeCompressed();
      
      }
      
      
      inline void print() const {
      
            std::cout << this << std::endl;
      
      }
      
      
      inline bool rowIsLocal(int globalRowIndex) const {
      
            return (this->global_starts_.first <= globalRowIndex && globalRowIndex < this->global_starts_.first + this->rows());
            
      }
      
      
      inline void setElement(_Scalar value, int localRowIndex, int localColumnIndex) {
      
            _Index r = (_Index) localRowIndex;
            _Index c = (_Index) localColumnIndex;
            if (this->coeff(r, c) != this->zero) {
                  this->coeffRef(r, c) = value;
            } else {
                  if (_Options == RowMajor) {
                        _Index newsizes [this->rows()];
                        for (_Index r = 0; r < this->row(); ++r) {
                              newsizes[r] = *(this->OuterIndexPtr() + r + 1) - *(this->OuterIndexPtr() + r);
                        }
                        newsizes[localRowIndex] += 1;
                        this->reserve(newsizes);
                        this->insert(r, c) = value;
                  } else {
                        _Index newsizes [this->cols()];
                        for (_Index c = 0; c < this->col(); ++c) {
                              newsizes[c] = *(this->OuterIndexPtr() + c + 1) - *(this->OuterIndexPtr() + c);
                        }
                        newsizes[localColumnIndex] += 1;
                        this->reserve(newsizes);
                        this->insert(r, c) = value;
                  }
                  this->makeCompressed();
            }
            
      }
      
      
      inline void setElementLocal(_Scalar value, int k) {
      
            _Index r = this->getLocalRowIndex(k);
            _Index c = this->getLocalColumnIndex(k);
            this->setElement(value, r, c);
      
      }
      
      
      inline void setZero(_Scalar zero) {
      
            this->zero = zero;
      
      }
      
      
      inline std::string toString() const {
      
            std::string result;
            result.append("Sparse Matrix A, id = " << this->id_ << ", rows = " << this->rows() << ", cols = " << this->cols() ", Storage = ");
            if (_Options == RowMajor) { 
                  result.append("Row-Major:" << std::endl);
            } else {
                  result.append("Col-Major:" << std::endl); 
            }
            for (int r = 0; r < this->rows(); ++r) {
                  for (_Index c = 0; c < this->cols(); ++c) {
                        result.append("A[" << r << "][" << c << "] = " << this->getElement(r, c) << std::endl);
                  }
            }            
            return result;
      
      }
      
      
      inline void initialize(int id, int proc, int nLocal, int mLocal, int i0, int j0, int rb, int cb) {
      
            this->setSubmatrixId(id);
            this->setProcId(proc);
            this->resize(nLocal, mLocal);
            this->global_starts_.first = i0;
            this->global_starts_.second = j0;
            this->setZero((_Scalar) 0.0);
            this->block_.first = rb;
            this->block_.second = cb;
      
      }
      
      
      inline void initialize(int id, int proc, int nLocal, int mLocal, int i0, int j0, int rb, int cb,
            const _Scalar * const * const matrix, bool copyGlobal) {
            
            this->setSubmatrixId(id);
            this->setProcId(proc);
            this->resize(nLocal, mLocal);
            this->global_starts_.first = i0;
            this->global_starts_.second = j0;
            this->setZero((_Scalar) 0.0);
            this->block_.first = rb;
            this->block_.second = cb;
            if (copyGlobal) {
                  for (_Index r = 0; r < this->rows(); ++r) {
                        for (_Index c = 0; c < this->cols(); ++c) {
                              this->(r, c) = matrix[r + i0][c + j0];
                        }
                  }
            } else {
                  for (_Index r = 0; r < this->rows(); ++r) {
                        for (_Index c = 0; c < this->cols(); ++c) {
                              this->(r, c) = matrix[r][c];
                        }
                  }
            }
            this->makeCompressed();
           
      }
      
      
      inline void initialize(int id, int proc, int nLocal, int mLocal, int i0, int j0, int rb, int cb,
            _Scalar value, int rowIndex, int colIndex) {
            
            this->setSubmatrixId(id);
            this->setProcId(proc);
            this->resize(nLocal, mLocal);
            this->global_starts_.first = i0;
            this->global_starts_.second = j0;
            this->setZero((_Scalar) 0.0);
            this->((_Index) rowIndex, (_Index) colIndex) = value;
            this->block_.first = rb;
            this->block_.second = cb;
            this->makeCompressed();
            
      }
      

